Implementation of an Automatic Slice-to-Slice Registration Tool

This paper outlines the basic steps in the design and implementation of a feature based Transmission Electron Microscopy (TEM) image registration application and highlights some of the implementation details, such as the detection of features, feature descriptor design, robust filtering of mismatched descriptors, and transform estimation. Although the approach chosen is based on the Scale Invariant Feature Transform (SIFT) method, it is optimized for the TEM image registration. Implementation of an automati sli e-to-sli e registration tool Pavel A. Koshevoy, Tolga Tasdizen, and Ross T. Whitaker April 27, 2006 Abstra t This paper outlines the basi steps in the design and implementation of a feature based Transmission Ele tron Mi ros opy (TEM) image registration appli ation and highlights some of the implementation details, su h as the dete tion of features, feature des riptor design, robust ltering of mismat hed des riptors, and transform estimation. Although the approa h hosen is based on the S ale Invariant Feature Transform (SIFT) method, it is optimizied for the TEM image registration. 1 Motivation The goal of this proje t is to provide a fully automati tool for sli e-to-sli e image registration of several hundred sli es assembled from high-resolution tile images. This tool is aimed at resear hers working with Transmission Ele tron Mi ros opy images. The hallenges lay in the fa t that ea h sli e is arbitrarily oriented in the imaging plane, and may have been warped independently from all other sli es. 2 Problem statement Given an ordered sequen e of sli es (e.g. S0, S1,... Sn) a transform must be onstru ted for ea h adja ent sli e pair that would map from the image spa e of sli e Si to the spa e of sli e Si+1. This task will be addressed within a feature mat hing framework. The problem an be partitioned into several sub-problems outlined below: • For ea h sli e, a gradient ve tor image pyramid and a Di eren e-of-Gaussian image pyramid must be onstru ted. • The extrema points of the DoG pyramid must be determined. • The dominant gradient ve tor orientation(s) in the neighborhood of ea h extrema point must be dete ted. • A des riptor for every dete ted gradient ve tor orientation of the extrema point must be generated. • For ea h pair of adja ent sli es, mat hing des riptors must be found. • Given the mat hing des riptors, a transform that best maps the extrema points from the image spa e of sli e A into the image spa e of sli e B must be al ulated. 3 Implementation details The spe i s of the onstru tion of the image pyramids are thoroughly overed by David G. Lowe[2℄ and will not be repeated here. Su e it to say, that a pyramid is a olle tion of o taves, where ea h o tave represents a redu tion of image resolution by a fa tor of 2. Ea h o tave is partitioned into a set of s ales where ea h su essive image is onvolved with a Gaussian lter of in reasing sigma value. The details of e ient implementation of this are overed by Lowe[2℄. 1 3.1 Dete ting extrema points The extrema points are the lo al minima and maxima points of the Di eren e-of-Gaussian image pyramids. Lowe[2℄ proposed looking for an extrema point in a 3×3×3 neighborhood within a DoG pyramid. However, experimentation has shown that this te hnique does not yield stri t extrema points that are greater than or less than all of the neighbors. Relaxing the extrema riteria to allow the extrema point to be equal to its neighbors yields a large number of adja ent extrema points. Therefore, an alternative method of extrema dete tion is proposed. Let D1 be a non-boundary image within the DoG pyramid. Let D0 be the image pre eeding D1 in the pyramid, and D2 the su eeding image. Assuming there are minima points within the D0, D1, D2 sli es of the pyramid, al ulate Amin = D0 − D1 Bmin = D2 − D1 The resulting images Amin and Bmin are thresholded to remap the negative values to zero. The minima point image is al ulated as Emin = Amin × Bmin The maxima point image is al ulated analogously. Amax = D1 − D0 Bmax = D1 − D2 Again, Amax and Bmax are thresholded to remap the negative values to zero. The maxima point image is al ulated as Emax = Amax × Bmax The resulting extrema point images Emin, Emax are thresholded to isolate strong maxima, and an 8onne ted lustering algorithms is used to dete t the peaks. For ea h luster, the key point is positioned at the enter-of-mass of the luster. 3.2 Dete ting des riptor orientations The des riptor has to be rotationally invariant, therefore it is ne essary to sele t a onsistent frame of referen e for sampling the neighborhood around the extrema point. The method that is urrently implemented in the appli ation follows the one des ribed by Lowe[2℄. Essentially, the neighborhood gradient orientation angles are a umulated into a 1D histogram. Ea h ontribution is weighed by the gradient magnitude and a 2D Gaussian weighting fun tion entered at the extrema point. The peaks of the histogram de ne the feature ve tor orientation angles. 3.3 Generating the des riptors During experimentation, several di erent des riptor generators were evaluated, in luding 2 versions of the des riptor re ommended by Lowe[2℄. All of them share the following properties: • The des riptors are based on extrema point neighborhood properties derived from the image (su h as the gradient ve tor image, or the extrema image). • The neighborhood is sampled within a lo al oordinate system based on the des riptor orientation angle. • The radius of the sampling window has to be large enough (in pixels) to apture the neighborhood properties. 2 The major di eren e between the alternate des riptor generators and the design proposed by Lowe rests in the way the sampling window is partitioned. Lowe re ommends that the des riptor onsist of a 4×4 ell grid of 8-bin gradient orientation histograms, whi h leads to a 128 dimensional des riptor ve tor. The downside of this design is that it dis ards information that falls outside the grid. The alternative design partitions the neighborhood into a set of on entri annuli, where ea h annulus is partitioned into a set of ells of equal area. Ea h ell may hold an orientation histogram as suggested by Lowe, or some other information (su h as the average extrema intensity values extra ted from Emin and Emax, or dominant gradient ve tor angle). Unfortunately, experimentation with alternative des riptor designs has not shown performan e improvement over the design proposed by Lowe. The performan e was evaluated in terms of the number of known mat hing des riptors being orre tly mat hed using brute for e mat hing. 3.4 Mat hing des riptors The mat hing pro ess is slightly di erent from the one outlined by Lowe. Lowe addresses a more general omputer vision problem, where dete tion of the same obje t at di erent s ales is important. The ele tron transmission mi ros opy images are typi ally taken at the same s ale, and undergo minor deformation on the global s ale, making the s ale invariant feature mat hing unne essary. Therefore, for the purposes of TEM image registration, the des riptors are mat hed against other des riptors sele ted from the same o tave and s ale of the pyramid. In order to a hieve s ale invarian e, all that is required is the mat hing of des riptors from any o tave and s ale of a pyramid against any other o tave and s ale of the other pyramid. This would, of ourse, in rease the number of mismat hes. A brute for e implementation of des riptor mat hing is not unreasonable for the purposes of this proje t. However, following in Lowes footsteps the urrent implementation uses an optimized kd-tree[5℄ with a bestbinrst nearest neighbor sear h algorithm[3℄. 3.5 Filtering out bad mat hes Lowe has suggested two ltering stages for removing poorly mat hed des riptors. The rst stage is based on the thresholding of the ratio of Eu lidian distan e (in des riptor spa e) between the query des riptor and its losest mat h to the distan e between the query des riptor and its se ond losest mat h. This is founded on the observation that a well mat hed des riptor is usually distin t enough from the se ond losest mat h that the ratio of distan es would fall below 0.5, where as the ratio of distan es for a mismat hed des riptor and its se ond losest mat h is typi ally greater than 0.5. Unfortunately, our experimental results on mat hing TEM images have shown that the ratio of the des riptor distan e ratio between losest and se ond losest mat h is not nearly as well separated for orre t mat hes and mismat hes, therefore this property an not be used for ltering out bad mat hes, as it dis ards pra ti ally all of the orre t mat hes as well. The se ond stage proposed by Lowe is based on lustering with the Hough transform[7℄, whi h will not be overed here. Su e it to say that in our implementation it was not as e e tive as the alternative method des ribed below. The performan e of the two lters was ompared in terms of the ratio of the dete ted orre t mat hes to the number of mat hes in the ltered set. An alternative lter that appears to be extremely e e tive for TEM images is based on the ratio of the distan e (in image spa e) between nearest extrema points in image Si, to the distan e between their mat hing points in the image Si+1. This lter relies on the assumption that the s ales of the images being mat hed are the same, whi h is true for the TEM images. Sin e the s ales are the same, the distan e between nearest neighbors in one image and the mat hing image should be nearly identi al. If the ratio of the two distan es deviates signi antly from 1.0, it an be assumed that one of the mat hes is wrong. When it is determined that one of the points is mismat hed, both of the mat hes are dis arded. The downside of this ltering approa h is that for every dis arded mismat h, it may also be dis arding a good mat h as well. 3.6 Estimating the transform The remaining set of mat hes may still ontain some mismat hes, whi h presents a problem for a Least Squares solution. Matthew Brown[4℄ proposed the use of RANSAC[6℄ to sele t a set of mat hes that de ne 3 a onsistent transform. Essentially, a few mat hes are sele ted at random to solve for the transform parameters. The number of initially sele ted mat hes depends on the number of transform parameters. For example, a 2nd order (a ne) bivariate Legendre polynomial transform has 6 parameters, it therefore requires 3 distin t mat hes. A 4th-order bivariate Legendre polynomial transform has 20 parameters, it requires 10 distin t mat hes. On e a transform has been estimated, the rest of the mat hes are veri ed as inliers or outliers. For ea h mat h point pair, the point expressed in the spa e of image Si is mapped via the transform into the spa e of image Si+1. The distan e of the mapped point to its mat h is used to lassify the mat h as an inlier or an outlier based on some threshold. The inliers and the original set of mat hes are then used to re-estimate the transform. This an be an iterative pro ess, where at ea h iteration the mat hes are lassi ed as inliers and outliers, until onvergen e or a maximum number of iterations is rea hed. Sin e the goal is to optimize the number of inliers, the pro ess is repeated with a new set of initial random mat hes, and the best results are kept. For further improvement, it is possible to sort the mat hes a ording to some metri , su h as the Eu lidean distan e between the des riptors in the des riptor spa e. Then, instead of uniform sampling, importan e sampling may be used to sele t initial mat hes for RANSAC. 3.7 Further re nement of the transform estimate Given a transform estimate, it may be possible to remat h the des riptors between the two images by restri ting the set of mat h andidates to a lo al neighborhood within the transform target image spa e. For example, an initial set of des riptor mat hes may be used to estimate a low order transform (e.g. a ne) between images Si and Si+1. Given the low order transform, ea h des riptor from image Si is mapped into image Si+1. Only the des riptors that fall within a lo al neighborhood of the mapped des riptor are onsidered for mat hing. This eliminates a number of potential mismat hes that would be in onsistent with the a ne transform. On e all the des riptors have been re-mat hed, RANSAC an be used again to estimate a higher order transform. 4 Results An example of typi al images that must be pro essed by our appli ation is given in gure 1 on the following page. A Di eren e-of-Gaussian and a gradient ve tor pyramid of 2 o taves with 3 s ales per o tave was onstru ted for ea h image. The extrema of the DoG pyramid are dete ted: 2951 points in the left image, 2953 points in the right image. For ea h dete ted extrema point the lo al gradient ve tor neighborhood is examined to determine dominant gradient ve tor orientations. For ea h dete ted orientation a des riptor is onstru ted. This results in 4732 des riptors in the left image, and 4601 des riptors in the right image. An illustration of the dete ted des riptors is given in gure 2 on the next page. The des riptors are mat hed resulting in 4601 mat hes. These mat hes are ltered down to 459 mat hes see gure 3 on page 6 for an illustration. RANSAC is used to sele t inliers onsistent with an a ne transform whi h results in 165 mat hes illustrated in gure 4 on page 6. The resulting registration is shown in gure 5 on page 7. Referen es [1℄ Lindeberg, T. 1994. S ale-spa e theory: A basi tool for analysing stru tures at di erent s ales. Journal of Applied Statisti s, 21(2):224-270. [2℄ Lowe, D.G. 2004. Distin tive Image Features from S ale-Invariant Keypoints. International Journal of Computer Vision. [3℄ Beis, J. and Lowe, D. G. 1997. Shape Indexing Using Approximate Nearest-Neighbour Sear h in HighDimensional Spa es. In Conferen e on Computer Vision and Pattern Re ognition, Puerto Ri o, pp. 10001006. 4 Figure 1: Sample sli es Two onse ution sli es. Ea h sli e was assembled from 12 high resolution Transmission Ele tron Mi ros opy images of a rabbit retina. Figure 2: The des riptors Visualization of the unmat hed des riptor ve tors dete ted in the two images: 4732 des riptors in the image on the left, 4601 on the right. 5 Figure 3: The ltered des riptor mat hes Visualization of the ltered mat hes 459 out of 4601 mat hes remain. Figure 4: RANSAC ltered mat hes Visualization of the RANSAC ltered mat hes 165 onsistent mat hes are sele ted out of 459 remaining mat hes. 6 Figure 5: sli e to sli e registration Visualization of the sli e-to-sli e registration results. 7 [4℄ Brown, M. and Lowe, D.G. 2002. Invariant Features from Interest Point Groups. In British Ma hine Vision Conferen e, Cardi , Wales, pp. 656-665. [5℄ Friedman, J.H., Bentley, J.L. and Finkel, R.A. 1977. An Algorithm for Finding Best Mat hes in Logarithmi Expe ted Time. ACM Transa tions on Mathemati al Software, 3(3):209-226. [6℄ Fis hler, M.A. and Bolles, R.C. 1981. Random Sample Consensus: A Paradigm for Model Fitting with Appli ations to Image Analysis and Automated Cartography. Communi ations of the ACM, 24(6):381395. [7℄ D. H. Ballard. 1981. Generalizing the Hough transform to dete t arbitrary patterns. Pattern Re ognition, 13(2):111-122.